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The role of solute attractive forces on hydrophobic interactions is studied by coordinated develop¬ 
ment of theory and simnlation results for Ar atoms in water. We present a concise derivation of the 
local molecnlar field (LMF) theory for the effects of solute attractive forces on hydrophobic interac¬ 
tions, a derivation that clarifies the close relation of LMF theory to the EXP approximation applied 
to this problem long ago. The simulation results show that change from purely repulsive atomic 
solute interactions to include realistic attractive interactions diminishes the strength of hydrophobic 
bonds. For the Ar-Ar rdfs considered pointwise, the numerical results for the effects of solute attrac¬ 
tive forces on hydrophobic interactions are of opposite sign and larger in magnitude than predicted 
by LMF theory. That comparison is discussed from the point of view of quasi-chemical theory, and it 
is suggested that the first reason for this difference is the incomplete evaluation within LMF theory 
of the hydration energy of the Ar pair. With a recent suggestion for the system-size extrapolation 
of the required correlation function integrals, the Ar-Ar rdfs permit evaluation of osmotic second 
virial coefficients B 2 - Those B 2 also show that incorporation of attractive interactions leads to more 
positive (repulsive) values. With attractive interactions in play, B 2 can change from positive to 
negative values with increasing temperatures. This is consistent with the historical work of Watan- 
abe, et al., that B 2 ~ 0 for intermediate cases. In all cases here, B 2 becomes more attractive with 
increasing temperature. 


I. INTRODUCTION 

The concept of a hydrophobic interaction is firmly em¬ 
bedded in general views of the folding of water-soluble 
protein molecules. Kauzmann [T] clearly articulated that 
idea: “Thus, proteins are stabilized by the same physical 
forces as those that keep oil and water from mixing ...” 
A key phenomenological point is that enthalpic hydration 
contributions to the thermodynamics of protein unfold¬ 
ing decrease, even vanish at moderate temperatures [T]. 
Hydrophobic interactions are then entropy dominated. 
The enthalpy-entropy balance depends importantly on 
temperature, and switches at higher temperatures [T]. 
The entropic hydration contributions to the thermody¬ 
namics of protein unfolding can vanish at higher temper¬ 
atures [T], and that condition has been called the “en¬ 
tropy convergence” point [5]. Nevertheless, below such an 
entropy convergence temperature, i.e., where hydropho¬ 
bic low-solubility is associated a negative entropy change, 
hydrophobic interactions get stronger with temperature 
increases though with reduced rate of strengthening. 
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From early days, the involvement of the hydration en¬ 
tropy has been conceptualized by imagining icebergs sur¬ 
rounding simple hydrophobic solutes, such as inert gases, 
e.g., Ar below. Tanford [3] attributed the original “ice¬ 
berg” language to G. N. Lewis. Silverstein [1], decades 
later, provides a modern view of the relevance of the ice¬ 
berg concept to hydrophobic phenomena. “Iceberg” is a 
widely recognized figure of speech, but has not been the 
basis of defensible statistical mechanical theory of these 
entropy effects. In fact, the statistical mechanical the¬ 
ory that eventually does explain the entropy convergence 
phenomena does not define or characterize iceberg struc¬ 
tures [2]- 

Because the iceberg parlance is vague and provocative, 
direct experimental demonstrations of so-called inverse- 
temperature behaviors are particularly helpful. Aggre¬ 
gation of sickle hemoglobin is a standard example [S]. 
Well-known aqueous polymers that separate with tem¬ 
perature increases, i.e., systems that exhibit a lower 
critical solution temperature (LOST), also provide ex¬ 
amples. Elastin-like peptides (ELPs) are probably the 
best known cases [MID]- Substantial molecular simula¬ 
tion work is available describing ELP collapse dlHIS] 
without addressing the statistical mechanical theory of 
hydrophobic interactions. Those descriptive simulation 
efforts are largely consistent with the traditional idea of 
the folding of elastin-like peptides upon heating, and with 
each other, but not entirely m- 
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Aqueous solutions of poly(N-isopropylacrylamide)s 
provide other examples of LCSTs [HI [20]. Polyethylene 
glycols (PEGs) in water also exhibit LCSTs [Hj. The 
polymers noted are water-soluble below their LCSTs. 
Thus they are substantially hydrophilic. But though they 
are complicated molecules, they directly demonstrate the 
inverse-temperature phenomena of classic hydrophobic 
effects. 

The successful statistical mechanical theory for the 
entropy convergence behavior developed over 

decades from counter-intuitive initial steps |5Dl - [55] . The 
statistical mechanical theory of hydrophobic interactions 
[35l [39] was formulated for hard sphere hydrophobic so¬ 
lutes in water, and theoretical progress has been as¬ 
sociated with attention to detail for such simple cases 
[241 EH 1^ . That methodical analysis strategy per¬ 
mits clarity in isolating the features that are the ul¬ 
timate interest. An important accomplishment of re¬ 
cent work [39] was then to prove numerically that rigor¬ 
ously defined hydrophobic interactions between atomic¬ 
sized hard sphere solutes in water also exhibit inverse- 
temperature behavior. Independently, new results for 
broader solute models arrived at consistent conclusions 
for those cases [40]. Building from those important ac¬ 
complishments, the present work investigates the theory 
for adding attractive inter-atomic forces to those primi¬ 
tive cases. 

The counter-intuitive ingredients of the statistical me¬ 
chanical theory [36] together with apparent disagreement 
with some experiments [SailTHls] that do involve attrac¬ 
tive forces, lead promptly to questions about the conse¬ 
quences of solute attractive forces associated with simple 
hydrophobic solutes [H]. That issue has been broadly 
discussed several times over the intervening years [4^144] - 
H5] without achieving a dehnitive solution. That situa¬ 
tion can now change on the basis of the new results for 
hydrophobic interactions noted above. 

Distinctions [49] deriving from inclusion of solute at¬ 
tractive forces are exemplified in FIG. and FIG.[^ In¬ 
clusion of solute attractive forces diminishes the strength 
of hydrophobic bonding: solvent attraction to the so¬ 
lute tends to pull the solute species apart. This behavior 
could be expected from sensitive appreciation [JU US] 
of preceding results. The local molecular field theory 
(LMF) discussed below is a simple, persuasive theory for 
these effects of attractive interactions m- Clarifying 
and testing that theory is the goal of this work. 

Though the substantially the same theory we test be¬ 
low was known and used [H] [ST] El] long ago, the in¬ 
sight underlying recent discussions of LMF theory, e.g., 
[50] . have considerably strengthened it. We here give a 
concise derivation with a clear analogy to a thermody¬ 
namic van der Waals picture, and is therefore unusually 
compelling. In the next section we outline the LMF the¬ 
ory. Numerical results, and conclusions are identified in 
Sec. |IV[ Methods for the several computational steps are 
collected in Sec. El 
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FIG. 1. Modeled radial distribution functions for WCA 
repulsive-force Ar solutes, based on the hard-sphere cavity 
distribution functions [39] . 
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FIG. 2. Ar-Ar radial distribution functions reconstructed 
from stratified (window) calculations. Notice (compare 
FIG. E that contact hydrophobic interactions are weaker 
when solute attractive forces are included. In contrast, 
solvent-separated correlations are more strongly structured 
with inclusion of atomic attractive forces. 


II. LOCAL MOLECULAR FIELD THEORY 

The LMF idea is to study the inhomogeneous density 
of a fluid subject to an external field. We focus on the 
density structure resulting from the placement of an Ar 
atom at a specific location. That distinguished atom ex¬ 
erts an external field on the surrounding fluid and distorts 
the density. With U the intermolecular potential energy 
function for the system and $ the external held exerted 
by the distinguished atom, the resulting distorted density 
is PaM (j’; U, $) at position r of a atoms of a molecule 
of type M. The goal of the LMF theory is to analyze 
PaM U, 4>) on the basis of the characteristics of the 
interactions U and <I>. 

We assume that a reference potential energy, de- 

























noted by has been identified to help in analyzing 

PaM (r; U, $). Specifically, our goal is the match 
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PaM (r;[/,$) = PaM (r; , (1) 

achieved for the reference system with interactions 
and an effective external field That effective field is 
the objective of the analysis below. A successful match 
Eq. § establishes aspects of U that can be treated as 
molecular mean-fields, thus offering a molecular mecha¬ 
nism for PaM (r', U, $). 

Identification of a reference potential energy function 
thus requires physical insight. One suggestion for 
the inter-atomic force fields of current simulation calcu¬ 
lations corresponds to Gaussian-truncated electrostatic 
interactions associated with the partial charges of simu¬ 
lation models m- In that case, the crucial difference 

U — (jraM — '>’7M'I) ) (2) 

aM,7M' 


is atom-pair decomposable. For the case of interest here, 

^aM 7 M' (l^’aM ” ’' 7 M' I) the WCA-attractive part of the 
Lennard-Jones interactions associated with the Ar atoms 

Ell¬ 
in seeking the match Eq. Q, we adopt an atom-based 
perspective, and focus on the chemical potential [S3], 


PaM = /3 Mn [paM {r; U, $) Aq,m^] 

+ PaM (r) + (r; p, f3U) , (3) 

of aM atoms, which decomposes 

$ = X! {rau) ■ (4) 

aM 

Here the temperature is T = (fcB/3)~^; the thermal de- 
Broglie wavelength A^m depends only on T and on fun¬ 
damental parameters associated with a atoms. As in¬ 
dicated, the excess contribution {r;p,j3U) depends 
functionally on {p,j3U), not on the external field. 

For some simulation models, the atom-based PaM 
(Eq. may raise questions regarding the operational 
status of atom chemical potentials. But this perspective 
would be satisfactory for ab initio descriptions of the so¬ 
lution, and is sufficiently basic that we do not further 
side-track this discussion. Similarly for the reference case 
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The forms Eqs. (|^ and ([^ allows us to express the match 
Eq. Q as 


</’aM (^) = ‘Po^M (r) 


PaM (r-, P, PU) - p^M^ (r; p, 


- rnnstpint 


( 7 \ 


The bracketed terms in Eq. Q depend functionally on 
the densities, not on the external field. The constant 
in Eq. 0 involves the chemical potentials of the two 
systems. 

The approximation 


PaM PU) ~ I^aM PU^°^) 

+ (r') w^^M'aM (W “ ^-1) dr' (8) 

7M' 

is then transparently analogous to van der Waals theory 
[54] for the inclusion of attractive interactions, i.e., Ap « 
—2ap, with a the van der Waals parameter describing 
attractive intermolecular interactions. Transcribing to 
the case of Ar(aq) at infinite dilution produces 


PAr (»’) ~ PAr (r) 

+ J [po (r') - Po] MoAr (I’’' - »’l) dr' • (9) 

The fields vanish far from their source, and therefore the 
constant contribution of Eq. Q is accommodated ex¬ 
plicitly in Eq. (|^. This argument matches the results 
of Rodgers and Weeks [S5j in the several instances they 
considered. Derivations that emphasize alternative (elec¬ 
trostatic) interactions are available elsewhere [5M55] . 

Though the statistical mechanical theory of Eq. 0 is 
simple, the field sought depends on the density, which 
depends on the field. A linear statistical mechanical ap¬ 
proximation Eq. (|^ produces the non-linear Eq. 0 to 
solve. The non-linearity is not an obstacle here because 
the densities on the right of Eq. ([^ are straightforwardly 
obtained from routine simulation (FIG.[^ see also [59] 1. 
Notice (FIG. that the effects of attractive ArO inter¬ 
actions on ArO correlations are modest, as was argued 
long ago [44] . 

Now consider pAr(j'; the density of Ar 

atoms without attractive interactions (r) but ex¬ 

periencing the effective field /Sp^r (’’) ■ We approximate 

m 


lnpAr(r';C/^°\$^°^)/pAr « “/^PAr (^) +ln?/Hs(»’) 

= -P {pAr (’’) - «aL (’’)) + In 5ArAr (^) > (10) 

adopting the repulsive-force solute results of FIG.[2 The 
field p^,) (r) incorporates aspects of the intermolecular 
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FIG. 3. Observed radial correlation of O atoms with an 
Ar atom, T = 300 K, p = 1 atm (heavy curve). Correlation 
functions (fainter, background curves) for hard-sphere model 
solutes with distances of closest approach 0.31 nm (FIG.[^ on 
the basis of cavity methods [39II59| . from Chaudhari [49]. The 
PC theory [35] predictions for the maxima of the hard sphere 
correlation functions would be close to 2, larger than these 
numerical results [SniEHl- For the soft-sphere case, attractive 
van der Waals interactions draw-in near-neighbor O-atoms 
slightly [311 EH . Since attractions draw-in, rather than draw- 
up, the attractive interactions case is not wetter than the 
reference case. 


ory for /lArO {r) is 


-In 


gArO jr) 

9Aroi'f'). 


hoo (r')po^MoAr(k 


-r|)dr', (14) 


where /iqo i^') is the observed 00 correlation function 
for pure water. Acknowledging closure approximations 
specific to traditional implementations, this is just the 
EXP approximation ISD applied to this correlation prob¬ 
lem long ago mi [S2I . This observation serves further to 


identify Eq. (13) as a relative of the EXP approximation 


also. Nevertheless, the distinction between the theory of 
Ref. I44| . with its specihc implementation details, from 
Eq. (13) should be kept in mind. The most prominent 


distinction is that Eq. (131 exploits liArO (’’) evaluated 


self-consistently or, here, the numerically exact result. 

Note further that the E q. ([I^ offers additional pos¬ 
sibilities compared to Eq. (]14[) for variety of outcomes 
because of possibilities from imbalance of (r) and 


-'"ArAr 


(r). 


B. Perspective from quasi-chemical theory 
14411481163] 


attractions as mean-field effects according to Eq. S- 
The match Eq. Q pairs this with 

InpArlr-;!/, ^j/pAr =lngArAr(»’) • (H) 


Combining with Eq. (101 


IngArAr (t-) = In ^ArAr (^) “ (f^Ar M (r)) 

- J [po {r') - po] dr' - r|) dr' . (12) 

Finally noting pAr = + ^ArAr ^nd rearranging 

yields 


— In 


gArAr (r) 

.gAiL (’’) 


/3«a\L (»’) 


+ J /iAro(r')po^MoAr(k'-»’l)dr'. (13) 


Quasi-chemical theory (QCT) provides insight into 
the LMF approximation Eq. ( jl^ . Since QCT is de¬ 
signed to evaluate interaction contributions to chemi¬ 
cal potentials,[3^inTl|M] Eq. Q is the relevant starting 
point. From the QCT formulation [48l[63], the packing 
contributions to those two chemical potentials are iden¬ 
tical, and cancel each other. Next to be considered [15] 
is the mean hydration energy, denoted by {e\r,n\ = 0), 
of the Ar appearing at r. That previous QCT effort [48] 
observed that the outer-shell QCT fluctuation contribu¬ 
tion was comparatively small. Thus (e|r, = 0) is the 

leading factor in describing the effect of attractive in¬ 
teractions being added [HI [48] . In the QCT study, [48] 
(e|r, 71^ = 0) was evaluated from molecular simulation 
data. The PC modeling of long-ago [U] recognized the 
importance of {e\r,n\ = 0), and used a RISM approxi¬ 
mation to incorporate the specific structure of the Ar 2 
diatom. Returning to the LMF theory, the right-most of 
Eq. Q addresses {e\r,n\ = 0), but does not calculate it 
for the detailed Ar 2 geometry. Incomplete evaluation of 
(e|r, nx = 0) is thus the chief neglect of the present LMF 
theory. 


A. Comparison to EXP theory |44| 


As noted above, the approximate theory Eq. (13) re¬ 
quires /lArO (^) J Slid we can conveniently take that from 
routine simulation. The corresponding approximate the- 


III. METHODS 
A. Simulations 

The simulations were carried-out with the GROMACS 
package [55], the SPC/E model of the water molecules 
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FIG. 4. Test of the LMF theory, Eq (13l. The net re¬ 
sult for the LMF approximation (black, dashed, right-side 
of Eq. (13l) is the sum of the direct interaction (blue, dot¬ 
ted curve) aud the water contribution (red, dotdashed curve, 

Eq. (^l). —In |^(?ArAr(r)/(;ArAr(’")] differs from the net LMF 


approximation both in contact and solvent-separated config¬ 
urations. 



r (nm) 


FIG. 5. Gomparison of LMF approximation Eq. (131 
with g^rAr (*') (reference system, fainter, dotted curve) and 


flArAr (r). Note the significantly different behavior of (JArAr (^) 
and qatAi (r) in the second shell, not addressed by this ap¬ 
proximation. 


[66] . and the OPLS force field. GROMACS selects SET¬ 
TLE constraint algorithm for rigid SPC/E water 
molecules. The same constraint algorithm was used in 
previous simulations involving water |391168] . Standard 
periodic boundary conditions were employed, with parti¬ 
cle mesh Ewald utilizing a cutoff of 1 nm and long-range 
dispersion corrections applied to energy and pressure. 
The Parrinello-Rahman barostat controlled the pressure 
at 1 atm, and the Nose-Hoover thermostat was used to 
maintain the temperature. The simulation cell for the 
Ar(aq) system consisted of two (2) argon molecules and 


1000 water molecules. Initial conhgurations were con¬ 
structed with PACKMOL [69] to construct a system close 
to the density of interest. The solute-solute separation 
spanning 0.33 nm to 1.23 nm was stratified using a stan¬ 
dard windowing approach and the results combined using 
the weighted histogram analysis method (WHAM) [70] . 
This involved 19 windows (and simulations) for window 
separations r ranging from 0.33 nm to 1.23 nm. 


B. Implementation of LMF theory 

With the information of FIG. [^ the LMF approxi¬ 
mation Eq. © depends only linearly on the attractive 
interactions. We evaluated Eq. © standardly, intro¬ 
ducing the spatial Fourier transforms 

*OAr W = J «OAr (^) ’ (1^) 

and 

hoAv (fc) = j hoAr (r) dr . (16) 

Then 

= J ^ArO W) PoPu^oAt (\r’ -r\) dr’. (17) 

The parameters for this application are po = 33.8/nm^, 
£OAr = 0.798 kJ/mol, coAr = 0.328 nm, eArAr = 
0.978 kJ/mol, and CTArAr = 0.340 nm. 


C. Osmotic i ?2 and Infinite Size Extrapolation 

The distribution function pArAr (’’) = ^ArAr (^) + l pro¬ 
vides access to the osmotic second virial coefficient, 

^2 = -^ lim //lArAr (r) d^r . (18) 

Z PAr-i-O J 

We utilize the extrapolation procedure of Kruger, et al. 

[num 

— 2 B 2 = lim 47r / /lArAr (r) tc(r/2i?)r^dr , (19) 

fl-i-oo Jq 

with 

w(x) = 1 - x-b . (20) 

Computed values for l/2i? > 0 were least-squares fitted 
with a polynomial quadratic order in 1 /2i?, then extrap¬ 
olated to l/2i? = 0. This procedure has been successfully 
tested [SHIITIIZI] and does not require further statistical 
mechanical theory. 
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FIG. 6. Extrapolation to evaluate the osmotic second virial 
coefficient, B 2 ^\ for the WCA repulsive-force Ar solutes 
(FIG. [^. T he sy mbol at l/2i? = 0 is the extrapolated 


value; see Sec III G Hydrophobic interactions gauged by 


become more attractive with increasing temperature in this 
range. 



75 

r 

1 


50 



/mol) 

25 

7 . 


CO 

s 


\ 

Vv 


0 


- ■■ 

cq' 





-25 


^ 


-50 

' 

\ 

_1_ 


0.5 


* 


* 


I 


* 280K 
V 300K 
■ 320K 

• 340K 
X 360K 


1.0 
1/2R (nm) 


1.5 


2.0 


FIG. 7. Extrapolation (Sec III G[ ) to evaluate the osmotic 
second virial coefficient, B 2 , for the full attractions case of 
FIG. Gomparing with FIG. we see that inclusion of 
solute attractive-forces makes these B 2 more positive (repul¬ 
sive). Hydrophobic interactions gauged by B 2 become more 
attractive with increasing temperature in this range. B 2 
changes from positive to negative values in T =(320K, 340K). 
Thus 52 ~ 0 in that interval, qualitatively consistent with the 
work of Watanabe, et al., 1451. 


IV. RESULTS AND DISCUSSION 

Changing purely repulsive atomic interactions to in¬ 
clude realistic attractions diminishes the strength of hy¬ 
drophobic bonds (FIGs. [^and[^. Within this LMF the¬ 
ory, and also the earliest theories for this [44l |48], the 
hydration environment competes with direct Ar-Ar at¬ 
tractive interactions (FIG.[^. The outcome of that com¬ 
petition is sensitive to the differing strengths of the at¬ 


tractive interactions. The earlier application |44j used 
the EXP approximation to analyze the available Monte 
Carlo calculations on atomic LJ solutes in water m 
That theoretical modelling found modest effects of at¬ 
tractive interactions, and encouraging comparison with 
the Monte Carlo results. This application of the LMF 
theory (Eq.[^ again predicts modest effects of attractive 
interactions, but the net comparison from the simulation 
results shows big differences. The ontcome alternative 
to the historical work is due to the fact that the earlier 
theory used the PC approximate results for the reference 
system g^rAr know that approximation 

is not accurate for this application [39] , despite being the 
only theory available. Here the LMF theory (Eq. (Hi) 
predicts modest-sized changes, though in addition it pre¬ 
dicts changes opposite in sign to the observed changes. 
Note further that g^rAr 9ArAr (j") differ distinc¬ 

tively in the second hydration shell, and those differences 
suggest more basic structural changes driven by attrac¬ 
tive interactions. 

The earliest study of these effects [31] went further to 
analyze a Lennard-Jones model of CH 4 -CH 4 (aq) with 
much larger attractive interactions. The theory devel¬ 
oped for that application was successful for the case stud¬ 
ied m. but the correspondence of that LJ model to CH 4 
solutes was not accurate enough |47j to warrant further 
interest. 

That earlier theory featured study of {e\r, n\ = 0) that 
has acquired a central role in QCT study of the present 
problem [48] . A more accurate evaluation would involve 
n-body (n > 2 ) correlations, perhaps even treated by nat¬ 
ural superposition approximations [76]. Detailed treat¬ 
ment of the Ar 2 diatom geometry is the most prominent 
difference between that QCT approach and the present 
LMF theory (Eq. @)- Nevertheless, a full QCT anal¬ 
ysis of these differences is clearly warranted and should 
be the subject of subsequent study. 

These changes due to attractive interatomic interac¬ 
tions are directly reflected in the values of B 2 (FIGs. 
and0. Slight curvature of the extrapolation (FIGs. 
andj^ is evident but, in view of the previous testing 
[6Hl|73l[74], not concerning. In all cases here, B 2 becomes 
more attractive with increasing temperature below T = 
360K. This behavior is consistent with accumulated expe¬ 
rience and recently obtained results (39] [dOj [73] . With at¬ 
tractive interactions in play, B 2 can change from positive 
to negative values with increasing temperatures. This is 
consistent with the historical work of Watanabe, et al., 
[45] that 52 « 0 for intermediate cases. 

Finally, we emphasize that since attractions make sub¬ 
stantial contributions, precise tests of the PG theory 
mm with results on cases with realistic attractive inter¬ 
actions should specifically address the role of attractive 
interactions that were not included in the PC theory. 
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